% MLE Estimation and Simulation

%need to specify
% run_estimate=0/1 
% run_simulate =0/1 simulate given current params
% run_opt_policy=0/1 solve for optimal policies and simulate

warning('off')


run_simulate=0
run_estimation=0
run_opt_policy=1
calc_std_error=0
run_subsidy_range=0
%ones below all need opt_policy
nat_level=0
cost_neut=0
min_damages=0
cost_half=0
new_store=0
no_reopt=0 % dont' recalculte optimal subs with new storage
net_cost=0
line_losses=0
trans_con=0 % transmission constraints
increase_renew=0 % takes 0, 1,2,3 for difference scenarios
clean_coal=0
marginal_increase=0 % marginal subsidy increases aroundq current system
tract_level=0   % optimal tract level subsidies. Can do min damages as well

over_dist_i=1 % households discount future subsidies at 15%
mcf125=0
mcf150=0

%don't turn on cost_neut if cost_half==1 or min_damages==1
%new_store needs opt_polciy (even though just gives gains of new_store
%given current subsidies
%to get optimal policy with new_store select new_store=1 cost_neut=1 (or
%min_damages=1

if min_damages==1
    cost_neut=0
end

if net_cost==1
    cost_neut=0
end


clearvars  -except run_estimation rev_neut run_simulate run_opt_policy calc_std_error ...
    run_subsidy_range nat_level cost_neut min_damages cost_half new_store net_cost line_losses increase_renew clean_coal ...
    marginal_increase tract_level no_reopt over_dist_i mcf125 mcf150 trans_con
close all; 
clc;

format long

% Global parameters-------------------------------------------------------

% ntract - Number of local labor markets.
% k - number of parameters estimated in first stage
% nstate - Number of states


global k ntract nstate nnerc nhours npp npol pol_types %nemis npol % nbins
global model_type  state_max statetaxcredit_t statemaxtaxcredit_t sunit_cap_t just_size net_metering C PV PV_hh line_losses_loc demand alpha1LL alpha2LL over_dist
global marginal_increase_global trans_con_loc
global scc


global tract_state_xxwalk_mat tract_nerc_xxwalk tract_state_xxwalk tract_nerc_xxwalk_mat %tract_nerc_xxwalk_mat tract_nerc_xxwalk_sparse
global  distance_border tract_border_xxwalk_mat population
  global  tl_global Sunit_tg Scost_tg Skwh_tg
  global PropTax IncTax SalesTax

    global PPIntCons PPRegionCons PPRegion2Cons PPInfCons PPRandStackedCons TransCap

%initialize parameters--------------------------------------------------
model_type='CollPercPolD' % NoDemD CollPercPolOwnD CollPercPolD CollPercPolIncD
%model_type='CollPercPolIncD'
scc=148; % 41 148   330


npol=4; % types of pollutants, co2 nox pm25 so2
net_metering=1; %one means some states don't have net metering and sale price differs from purchase price
C=.02;
over_dist=over_dist_i;
%k=1 sigma
%k=2 mean_gamma
%k=next =demographics
%k=end-1, end, gamma_N and gamma_2n
NCopies=2;

pol_types=["co2"; "nox"; "pm25"; "so2"];

marginal_increase_global=0;
tl_global=0;
Sunit_tg=0;
Scost_tg=0;
Skwh_tg=0;

if strcmp(model_type,'NoDem')==1 || strcmp(model_type,'NoDemD')==1
    k=4; % # parameters in utiltiy function
elseif strcmp(model_type,'CollPerc')==1 || strcmp(model_type,'CollPercD')==1 % college percentage in utiilty
        k=5;
elseif strcmp(model_type,'CollPercPol')==1 || strcmp(model_type,'CollPercPolD')==1 % college and percent Dem
        k=6 ;  
elseif strcmp(model_type,'CollPercPolOwn')==1 || strcmp(model_type,'CollPercPolOwnD')==1 || strcmp(model_type,'CollPercPolIncD')==1 % college and percent Dem and percent owner occupied
        k=7 ;          
end



if over_dist==1
    model_type=strcat(model_type, "OD");
%    spec=strcat(spec, "CC");
end

PV=19.39792919; % 25 years with 2% interest rate and .5% depreciation

%PV_hh= 9.704591; % 25 years with 10% interest rate and .5% depreciation

PV_hh= 7.247318; % 25 years with 15% interest rate and .5% depreciation

nhours=24*365-10;

% Read in data-----------------------------------------------------------



[A_t, M_t, BI_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, tract_nerc_xxwalk, tract_state_xxwalk,  distance_border, tract_border_xxwalk, population]=read_tract_data3();

ntract=size(A_t,1);


tract_nerc_xxwalk_sparse=sparse(tract_nerc_xxwalk,1:size(tract_nerc_xxwalk,1),1)';
tract_nerc_xxwalk_mat=full(tract_nerc_xxwalk_sparse);

tract_border_xxwalk_sparse=sparse(tract_border_xxwalk,1:size(tract_border_xxwalk,1),1)';
tract_border_xxwalk_mat=full(tract_border_xxwalk_sparse);


[Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state, PIns0_state, statemaxtaxcredit, statetaxcredit, sunit_cap, StateNames]=...
    read_state_data();

nstate=size(Sunit_state,1);



 title='../Data/CleanedData/CountySunlightProfileT.csv';
CountySunlightProfile=readmatrix(title);  

CountySunlightProfile=CountySunlightProfile';

CountySunlightProfile=CountySunlightProfile./(sum(CountySunlightProfile,2)*ones(1,nhours)); % normalize to sum to one. over course of year. 
%different than write up but easier for coding

TractSunlightProfile=zeros(ntract,nhours);

 title='../Data/CleanedData/county_i.csv'; % owner occupied perc. Normalized to mean 0
county_i_xxwalk=readmatrix(title);


for t_i=1:ntract
    c_i=county_i_xxwalk(t_i,1);
    TractSunlightProfile(t_i, :)=CountySunlightProfile(c_i,:);
end


tract_state_xxwalk_mat=sparse(tract_state_xxwalk,1:size(tract_state_xxwalk,1),1)';

tract_state_xxwalk_mat=full(tract_state_xxwalk_mat); %ntract by nstate

nnerc=max(tract_nerc_xxwalk);




% convert state data to tract data

Sunit_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit_state'),2);
Scost_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost_state'),2);
Skwh_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh_state'),2);

nbin=size(Nbar_i,2);

statemaxtaxcredit_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*statemaxtaxcredit'),2)*ones(1,nbin);
statetaxcredit_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*statetaxcredit'),2)*ones(1,nbin);
sunit_cap_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*sunit_cap'),2)*ones(1,nbin);

P_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*P_state'),2);
PSale_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*P_sale_state'),2);
PIns_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*PIns_state'),2);
PIns0_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*PIns0_state'),2);

[LoadRaw, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf, PPScale, ...
  ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2]=read_elec_data();


rng(123)


PPRandStacked=zeros(npp,nhours,NCopies);

for i=1:NCopies
    PPRandStacked(:,:,i)=(PPScale*ones(1,nhours)).*randn(npp,nhours);
end





ResProd=calc_res_prod(A_t, M_t,TractSunlightProfile); % this uses data level installations. Could also do with predicted installations


Load=LoadRaw+ResProd; % load raw doesn't account for residential 

ELoadData=LoadRaw-RenProd;

ELoadMax=max(ELoadData,[],1);

ELoadMax99=prctile(ELoadData,99,1);

line_losses_loc=0;
if line_losses==1

    line_losses_loc=1;

    title='../Data/electricity-generation/line-loss-alphas.csv';

    alphas=readmatrix(title);

    alpha1LL=alphas(:,2);

    alpha2LL=alphas(:,3);


    %solve for raw demand using quadratic equation
    a_temp=Load-alpha1LL';
    b_temp=alpha2LL';
    c_temp=ResProd;


    demand=sqrt(4.*a_temp.*b_temp-4.*b_temp.*c_temp+1)+2.*b_temp.*c_temp-1;

    demand=demand./(2.*b_temp);

    Load_test=demand+calc_line_losses(demand,ResProd);

end



trans_con_loc=0;
if trans_con==1

    trans_con_loc=1;

    title='../Data/CleanedData/PPInt-transmission-low.csv';
    PPInt=readmatrix(title);

    title='../Data/CleanedData/PPScale-transmission-low.csv';
    PPScale=readmatrix(title);

    PPRegion=zeros(npp, nnerc);

    PPRegion2=zeros(npp,nnerc);



PPInf=zeros(npp,nnerc);


for r_i=1:nnerc
   
    title=(sprintf('../Data/CleanedData/PPRegion-transmission-low_r%d.csv',r_i));    
%    title=(sprintf('../Data/CleanedData/PPRegion-own-region_r%d.csv',r_i));    
    PPRegion(:, r_i)=readmatrix(title);   
    
    title=(sprintf('../Data/CleanedData/PPRegion2-transmission-low_r%d.csv',r_i));  
 %   title=(sprintf('../Data/CleanedData/PPRegion2-own-region_r%d.csv',r_i));  
    PPRegion2(:, r_i)=readmatrix(title);     
    
    title=(sprintf('../Data/CleanedData/PPInf-transmission-low_r%d.csv',r_i));  
%    title=(sprintf('../Data/CleanedData/PPInf-own-region_r%d.csv',r_i));  
    PPInf(:, r_i)=readmatrix(title);         
    
end




     title='../Data/CleanedData/PPInt-transmission-high.csv';
    PPIntCons=readmatrix(title);
    

    title='../Data/CleanedData/PPScale-transmission-high.csv';
    PPScaleCons=readmatrix(title);


    PPRegionCons=zeros(npp, nnerc);

    PPRegion2Cons=zeros(npp,nnerc);
    PPInfCons=zeros(npp,nnerc);

    TransCap=zeros(nnerc,1);

    for r_i=1:nnerc
    %    title=(sprintf('../Data/CleanedData/PPRegionCons_r%d.csv',r_i));    
        title=(sprintf('../Data/CleanedData/PPRegion-transmission-high_r%d.csv',r_i));    
        PPRegionCons(:, r_i)=readmatrix(title);   
        
    %    title=(sprintf('../Data/CleanedData/PPRegion2Cons_r%d.csv',r_i));  
        title=(sprintf('../Data/CleanedData/PPRegion2-transmission-high_r%d.csv',r_i));  
        PPRegion2Cons(:, r_i)=readmatrix(title);   
        
    
         title=(sprintf('../Data/CleanedData/PPInf-transmission-high_r%d.csv',r_i));  
         PPInfCons(:, r_i)=readmatrix(title);      


         TransCap(r_i,1)=prctile(ELoadData(:,r_i), 75);
    
    end

    PPRandStackedCons=zeros(npp,nhours,NCopies);
    
    for i=1:NCopies
        PPRandStackedCons(:,:,i)=(PPScaleCons*ones(1,nhours)).*randn(npp,nhours);
    end


end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%Begin Estimation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Estimate model----------------------------------------------------------
if run_estimation==1 
    just_size=0;  %turn on to only target moments regarding size of array
    %%%% investigate correlations before estimation to get starting values
    Prob_b_data=BI_t./ sum(HH_bin,2);
    
    HH_t= sum(HH_bin,2);

    N_mean=15;

    X=((P_t+Skwh_t).*A_t.*N_mean-(1-Scost_t).*(PIns0_t+ PIns_t.*N_mean)+N_mean.*Sunit_t);
    

    %test3=regress(log(Prob_b_data(Prob_b_data>0)), [X(Prob_b_data>0) CollPerc_t(Prob_b_data>0) Pol_t(Prob_b_data>0) ones(size(P_t(Prob_b_data>0))) ] );
    guess_N=7.5;
    guess_2N=-0.25;
    
    if strcmp(model_type,'CollPercPol')==1 || strcmp(model_type,'CollPercPolD')==1 || strcmp(model_type,'CollPercPolDOD')==1
        test4=lscov([X(Prob_b_data>0) CollPerc_t(Prob_b_data>0) Pol_t(Prob_b_data>0) ones(size(P_t(Prob_b_data>0))) ],...
            log(Prob_b_data(Prob_b_data>0)), (HH_t(Prob_b_data>0)) );

        guess_sigma=1./test4(1,1)
        guess_col=test4(2,1).*guess_sigma;
        guess_pol=test4(3,1).*guess_sigma;

    %      params_guess=[guess_sigma; -100; guess_col; guess_pol; guess_N; guess_2N]    
          params_guess=[guess_sigma; -100; guess_col; guess_pol; guess_N; guess_2N]    
    elseif  strcmp(model_type,'CollPercPolOwn')==1 || strcmp(model_type,'CollPercPolOwnD')==1
        test4=lscov([X(Prob_b_data>0) CollPerc_t(Prob_b_data>0) Pol_t(Prob_b_data>0) Own_t(Prob_b_data>0) ones(size(P_t(Prob_b_data>0))) ],...
        log(Prob_b_data(Prob_b_data>0)), (HH_t(Prob_b_data>0)) );

        guess_sigma=1./test4(1,1)
        guess_col=test4(2,1).*guess_sigma;
        guess_pol=test4(3,1).*guess_sigma;
        guess_own=test4(4,1).*guess_sigma;

          params_guess=[guess_sigma; -100; guess_col; guess_pol; guess_own; guess_N; guess_2N] 

    elseif  strcmp(model_type,'CollPercPolInc')==1 || strcmp(model_type,'CollPercPolIncD')==1
        test4=lscov([X(Prob_b_data>0) CollPerc_t(Prob_b_data>0) Pol_t(Prob_b_data>0) Inc_t(Prob_b_data>0) ones(size(P_t(Prob_b_data>0))) ],...
        log(Prob_b_data(Prob_b_data>0)), (HH_t(Prob_b_data>0)) );

        guess_sigma=1./test4(1,1)
        guess_col=test4(2,1).*guess_sigma;
        guess_pol=test4(3,1).*guess_sigma;
        guess_inc=test4(4,1).*guess_sigma;

          params_guess=[guess_sigma; -100; guess_col; guess_pol; guess_inc; guess_N; guess_2N] 
     elseif  strcmp(model_type,'CollPerc')==1 || strcmp(model_type,'CollPercD')==1
        test4=lscov([X(Prob_b_data>0) CollPerc_t(Prob_b_data>0)  ones(size(P_t(Prob_b_data>0))) ],...
        log(Prob_b_data(Prob_b_data>0)), (HH_t(Prob_b_data>0)) );

        guess_sigma=1./test4(1,1)
        guess_col=test4(2,1).*guess_sigma;


          params_guess=[guess_sigma; -100; guess_col; guess_N; guess_2N] 
      elseif  strcmp(model_type,'NoDem')==1 || strcmp(model_type,'NoDemD')==1
        test4=lscov([X(Prob_b_data>0)   ones(size(P_t(Prob_b_data>0))) ],...
        log(Prob_b_data(Prob_b_data>0)), (HH_t(Prob_b_data>0)) );

        guess_sigma=1./test4(1,1)
       % guess_col=test4(2,1).*guess_sigma;


          params_guess=[guess_sigma; -100;  guess_N; guess_2N] 
    end

%title=(sprintf('Outputs/params_est%snm%d.csv',model_type, net_metering));
%  params_guess=readmatrix(title);

    state_max=1; % subsidies are subject to state level maximums
    
    gmm_dist=calc_gmm3(params_guess, A_t,  M_t, BI_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
        Sunit_t, Scost_t, Skwh_t, P_t, PSale_t, PIns_t, PIns0_t)
    
    
    options = optimset...
        ('Display','iter-detailed', 'MaxIter', 1000, 'MaxFunEvals', 1000, 'TolFun', 1E-7);
    
      


%% iterate over selected params
 

 params_change=[2;3;4;5;6]; % params I want to iterate over. ALl others are fixed. 

 params_change=ones(k-1,1);

for i=1:length(params_change)
    params_change(i,1)=i+1; 
end 

param_locs=zeros(k,1); 
for i=1:length(params_change)
    param_locs(params_change(i),1)=1; 
end

params_est1=params_guess;

    f = @(params_frag)calc_gmm3...
        (params_est1.*(1-param_locs) + param_locs.*[0; params_frag], A_t,  M_t, BI_t, Nbar_i,  HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
     Sunit_t, Scost_t, Skwh_t, P_t, PSale_t,PIns_t, PIns0_t);

    [params_frag_est, fval, exitflag, output]...
        = fminunc(f,params_est1(params_change,1),options);
    
    params_est15=params_est1.*(1-param_locs) + param_locs.*[0; params_frag_est]
    % constant   
    %}

    
        gmm_dist=calc_gmm3(params_est15, A_t,  M_t, BI_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
        Sunit_t, Scost_t, Skwh_t, P_t,PSale_t, PIns_t, PIns0_t)
    
    
%% all params

    f = @(params)calc_gmm3...
        (params, A_t,  M_t, BI_t, Nbar_i,  HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t,...
     Sunit_t, Scost_t, Skwh_t, P_t, PSale_t,PIns_t, PIns0_t);

    [params_est, fval, exitflag, output]...
        = fminunc(f,params_est15,options);

    params_est
    
    [params_est, fval, exitflag, output]...
        = fminsearch(f,params_est,options); 
    

    params_est
    
    gmm_dist=calc_gmm3(params_est, A_t,  M_t, BI_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
        Sunit_t, Scost_t, Skwh_t, P_t,PSale_t, PIns_t, PIns0_t)
    
    
  

   title=(sprintf('Outputs/params_est%snm%d.csv',model_type, net_metering));
%      csvwrite(title, params_est);
       writematrix(params_est, title);

end


    if calc_std_error==1

        %sloppy: a few variables passed as globals, not funciton arguemnts:
        %global distance_border tract_border_xxwalk_mat population
        %global PropTax IncTax SalesTax

        PropTax_o=PropTax;
        IncTax_o=IncTax;
        SalesTax_o=SalesTax;
        distance_border_o=distance_border;
        tract_border_xxwalk_mat_o=tract_border_xxwalk_mat;
        population_o=population;

title=(sprintf('Outputs/params_est%snm%d.csv',model_type, net_metering));
  params_est=readmatrix(title);

        n_boot=100;


            options = optimset...
        ('Display','final', 'MaxIter', 1000, 'MaxFunEvals', 1000, 'TolFun', 1E-7);

        length=size(A_t,1);

        params_boot=zeros(k,n_boot);

        for boot=1:n_boot
            order=randi(length, length,1);


             PropTax=PropTax_o(order);
             IncTax=IncTax_o(order);
             SalesTax=SalesTax_o(order);
             distance_border=distance_border_o(order);
             tract_border_xxwalk_mat=tract_border_xxwalk_mat_o(order,:);
             population=population_o(order);

                f = @(params)calc_gmm3...
            (params, A_t(order),  M_t(order), BI_t(order), Nbar_i(order,:),  HH_bin(order,:), CollPerc_t(order), Pol_t(order), Own_t(order), Inc_t(order), ...
         Sunit_t(order), Scost_t(order), Skwh_t(order), P_t(order), PSale_t(order), PIns_t(order), PIns0_t(order));

        [params_boot(:,boot), fval, exitflag, output]...
            = fminsearch(f,params_est,options);

             PropTax=PropTax_o;
             IncTax=IncTax_o;
             SalesTax=SalesTax_o;
             distance_border=distance_border_o;
             tract_border_xxwalk_mat=tract_border_xxwalk_mat_o;
             population=population_o;

        end

        std_error=std(params_boot,0,2);

       title=(sprintf('Outputs/std_error%snm%d.csv',model_type, net_metering));
         % csvwrite(title, std_error);  
          writematrix(std_error, title );  

    end
 
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%Begin Simulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Read in parameters-------------------------------------------------------

title=(sprintf('Outputs/params_est%snm%d.csv',model_type, net_metering));
  params=readmatrix(title);



%% Baseline Simulation 
if run_simulate==1
    
    state_max=1 % subject to state level maximum subsidies from data
    
    run_damages=1;
    



    [M_sim,BI_sim, D_sim, DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);




     ResProd=calc_res_prod(A_t, M_sim, TractSunlightProfile); %yearly res prod in mwh

     marg_damages_t=calc_marg_damages4(MargDamagesHour, TractSunlightProfile, ResProd); % this is per year. multiplied by A_t later to get lifetime



%total damages avoided per year by tract (compare to Sexton
%et al, roughly 200 to 1200 dollars. Figure 3 1225). multiply by 16 because
%theirs is 4w system ~ 16 panels
TDA=16.*A_t.*marg_damages_t./PV;

region_ave=zeros(nnerc,1);

for n_i=1:nnerc
    region_ave(n_i,1)=-1000*sum(TDA.*(tract_nerc_xxwalk==n_i))./sum(tract_nerc_xxwalk==n_i);
    region_A_t(n_i,1)=16*sum(A_t.*(tract_nerc_xxwalk==n_i))./sum(tract_nerc_xxwalk==n_i)./PV;
end

state_ave=zeros(nstate,1);

for s_i=1:nstate
    state_ave(s_i,1)=-1000*sum(TDA.*(tract_state_xxwalk==s_i))./sum(tract_state_xxwalk==s_i);
end


state_perkwh=zeros(nstate,1);
region_ave


for s_i=1:nstate
    state_perkwh(s_i,1)=-1000*sum(marg_damages_t.*(tract_state_xxwalk==s_i))./sum(tract_state_xxwalk==s_i);
end

nnerc_perkwh=zeros(nnerc,1);

for n_i=1:nnerc
    nnerc_perkwh(n_i,1)=-1000*sum(marg_damages_t.*(tract_nerc_xxwalk==n_i))./sum(tract_nerc_xxwalk==n_i);
end


    Prob_b_sim=BI_sim./ sum(HH_bin,2);   
    
  %  total_panels_base=sum(M_sim);
   total_installs_base=sum(BI_sim);
    
    title=(sprintf('Outputs/M_t%snm%d%s.csv', model_type, net_metering, "Base"));
    %csvwrite(title, M_sim); 
    writematrix(M_sim, title);

    

    title=(sprintf('Outputs/D%snm%d%s.csv', model_type, net_metering, "Base")); % damages
   % csvwrite(title, D_sim); 
    writematrix(D_sim, title);
    
    title=(sprintf('Outputs/BI%snm%d%s.csv', model_type, net_metering, "Base")); % buildings
  %  csvwrite(title, BI_sim);  
  writematrix(BI_sim, title);
    
    title=(sprintf('Outputs/Fisc%snm%d%s.csv', model_type, net_metering, "Base")); % buildings
    %csvwrite(title, Exp_sub_sim);      
    writematrix(Exp_sub_sim, title);

    title=('Outputs/HH_tract.csv'); % households (google)
  %  csvwrite(title, sum(HH_bin,2));    
  writematrix(sum(HH_bin,2), title);

  %once to get damages offset (so damages with no rooftop solar
[D0, DPol0, ~]=calc_damages5(0*ResProd, Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ProdMedian, ...
    DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2); % damages\\


    title=(sprintf('Outputs/D0%s.csv', model_type  )); % damages
   % csvwrite(title, D_sim); 
    writematrix(D0, title);


%Crago Chernyakhovskiy
    CCstates=[6; 7; 17; 18; 19; 27; 28; 30; 36; 37; 43; 46];
    Sunit_stateCC=Sunit_state;
    Sunit_stateCC(CCstates)=Sunit_state(CCstates)+.25; %250 watts per panel but all is in 1000D
    
%     [M_simCC,BI_simCC, D_simCC,M_bin_simCC, pi_binCC, Exp_sub_simCC]=simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, ...
%    Sunit_stateCC, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state,PIns0_state); 


    [M_simCC,BI_simCC, D_simCC,  DPol_simCC, MargDamagesHourCC, M_bin_simCC, pi_binCC, Exp_sub_simCC,  Tot_ut_simCC]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_stateCC, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

    BI_sim_state=sum((BI_sim*ones(1,nstate)).*tract_state_xxwalk_mat)';  
    BI_sim_stateCC=sum((BI_simCC*ones(1,nstate)).*tract_state_xxwalk_mat)'; 
    
    BI_sim_tot=sum(BI_sim_state(CCstates));
    BI_sim_totCC=sum(BI_sim_stateCC(CCstates));

CC_change=(BI_sim_totCC-BI_sim_tot)./BI_sim_tot

%gillinagham tsvetanov

eps=.00001;

    [M_simGT,BI_simGT, D_simGT,  DPol_simGT, MargDamagesHourGT, M_bin_simGT, pi_binGT, Exp_sub_simGT,  Tot_ut_simGT]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state,PIns0_state+eps, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

  

    elast_global=(sum(BI_simGT-BI_sim)./sum(BI_sim)).*(20./eps)



    stop
 
% model-based decomposition
%1: equalize subsidies


state_max=0; % turn off state maximums because they differ across states

net_metering_o=net_metering; 
net_metering=0; % assume all states offer net metering

HH_tract= sum(HH_bin,2);

HH_state=(HH_tract'*tract_state_xxwalk_mat)';

Sunit_state_cons=ones(size(Sunit_state))*(sum(Sunit_state.*HH_state)./sum(HH_state));
Skwh_state_cons=ones(size(Sunit_state))*(sum(Skwh_state.*HH_state)./sum(HH_state));
Scost_state_cons=ones(size(Sunit_state))*(sum(Scost_state.*HH_state)./sum(HH_state));

mult=1; % scale by certain value to make total installations equal to baseline
position=1; %which object to scale? 1 is subsidies

%[distance]=calc_decomp_dist(mult, position, total_installs_base,params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, ...
%Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state, P_sale_state, PIns_state,PIns0_state); 

options = optimset...
      ('Display','iter-detailed', 'MaxIter', 30, 'MaxFunEvals', 30, 'TolFun', 1E-3);
    
    f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

    [mult, fval, exitflag, output]...
        = fminunc(f,mult,options);
    
Sunit_state_cons=Sunit_state_cons.*mult;
Skwh_state_cons=Skwh_state_cons.*mult;
Scost_state_cons=Scost_state_cons.*mult;  

%[M_sim,BI_sim, D_sim,M_bin_sim, pi_bin, Exp_sub_sim]=simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, ...
%Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state, P_sale_state, PIns_state,PIns0_state); % will also need damages stuff

    [M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);


 Prob_b_NS=BI_sim./ sum(HH_bin,2);  
 


title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NS"));
%csvwrite(title, M_sim); 
writematrix(M_sim, title)
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NS")); % damages
writematrix(D_sim, title) 
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NS")); % buildings
writematrix(BI_sim, title) 
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NS")); % Fiscal cost
%csvwrite(title, Exp_sub_sim);
writematrix(Exp_sub_sim, title)





%2. Equalize electricity prices

P_state_cons=ones(size(P_state))*(sum(P_state.*HH_state)./sum(HH_state));
P_sale_state_cons=ones(size(P_state))*(sum(P_sale_state.*HH_state)./sum(HH_state));

mult=1; % scale by certain value to make total installations equal to baseline
position=2; %which object to scale? 2 is elec prices

 f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

    [mult, fval, exitflag, output]...
        = fminunc(f,mult,options);
    
P_state_cons=P_state_cons.*mult;
P_sale_state_cons=P_sale_state_cons.*mult;
 
    
%[M_sim,BI_sim, D_sim,M_bin_sim, pi_bin, Exp_sub_sim]=simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, ...
%Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state,PIns0_state); % will also need damages stuff


    [M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NSP"));
%csvwrite(title, M_sim); 
 writematrix(M_sim, title);
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NSP")); % damages
%csvwrite(title, D_sim); 
 writematrix(D_sim, title);
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NSP")); % buildings
 writematrix(BI_sim, title);
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NSP")); % Fiscal cost
%csvwrite(title, Exp_sub_sim);
 writematrix(Exp_sub_sim, title);

%3. Equalize installation prices
PIns_state_cons=ones(size(PIns_state))*(sum(PIns_state.*HH_state)./sum(HH_state));
PIns0_state_cons=ones(size(PIns0_state))*(sum(PIns0_state.*HH_state)./sum(HH_state));

mult=1; % scale by certain value to make total installations equal to baseline
position=3; %which object to scale? 2 is elec prices

  
 f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state_cons,PIns0_state_cons, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

    [mult, fval, exitflag, output]...
        = fminunc(f,mult,options);
    
PIns_state_cons=PIns_state_cons.*mult;
PIns0_state_cons=PIns0_state_cons.*mult;

    [M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state_cons,PIns0_state_cons, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NSPI"));
 writematrix(M_sim, title);
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NSPI")); % damages
 writematrix(D_sim, title); 
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NSPI")); % buildings
 writematrix(BI_sim, title);  
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NSPI")); % Fiscal cost
 writematrix(Exp_sub_sim, title);

%4. Equalize sunslight
A_t_cons=ones(size(A_t))*(sum(A_t.*HH_tract)./sum(HH_tract));


mult=1; % scale by certain value to make total installations equal to baseline
position=4; %which object to scale? 4 is sunlight
 
 f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t_cons, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state_cons,PIns0_state_cons, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

[mult, fval, exitflag, output]...
    = fminunc(f,mult,options);
    
A_t_cons=A_t_cons.*mult;


[M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
    simulate(params, A_t_cons, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state_cons,PIns0_state_cons, ...
Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);


title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NSPIS"));
 writematrix(M_sim, title);
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NSPIS")); % damages
 writematrix(D_sim, title); 
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NSPIS")); % buildings
 writematrix(BI_sim, title); 
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NSPIS")); % Fiscal cost
 writematrix(Exp_sub_sim, title);

%5. Equalize dems
CollPerc_t_cons=ones(size(A_t))*(sum(CollPerc_t.*HH_tract)./sum(HH_tract));
Pol_t_cons=ones(size(A_t))*(sum(Pol_t.*HH_tract)./sum(HH_tract));
Own_t_cons=ones(size(A_t))*(sum(Own_t.*HH_tract)./sum(HH_tract));
Inc_t_cons=ones(size(A_t))*(sum(Inc_t.*HH_tract)./sum(HH_tract));

%equalize bin sizes
%{
HH_bin_dist=sum(HH_bin)./sum(sum(HH_bin)); % mean proportion in each bin
HH_bin_eq=HH_tract.*HH_bin_dist;

Nbar_eq=Nbar_i(1,:).*ones(ntract,1);
%}

mult=1; % scale by certain value to make total installations equal to baseline
position=5; %which object to scale? 5 is dem

 f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t_cons, Nbar_i, HH_bin, CollPerc_t_cons, Pol_t_cons, Own_t_cons, Inc_t_cons, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state_cons,PIns0_state_cons, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

[mult, fval, exitflag, output]...
    = fminunc(f,mult,options);
    
CollPerc_t_cons=CollPerc_t_cons.*mult;
Pol_t_cons=Pol_t_cons.*mult;
Own_t_cons=Own_t_cons.*mult;

[M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
    simulate(params, A_t_cons, Nbar_i, HH_bin, CollPerc_t_cons, Pol_t_cons, Own_t_cons, Inc_t_cons, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state_cons,PIns0_state_cons, ...
Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NSPISD"));
 writematrix(M_sim, title); 
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NSPISD")); % damages
 writematrix(D_sim, title);
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NSPISD")); % buildings
 writematrix(BI_sim, title);
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NSPISD")); % Fiscal cost
 writematrix(Exp_sub_sim, title);  


%5 XXXX do electricty prices only
net_metering=net_metering_o; %reset to what it was before for decomp

mult=1; % scale by certain value to make total installations equal to baseline
position=2; %which object to scale? 2 is elec prices

 f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state, Skwh_state, Scost_state, P_state_cons, P_sale_state_cons, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

    [mult, fval, exitflag, output]...
        = fminunc(f,mult,options);
    
P_state_cons=P_state_cons.*mult;
P_sale_state_cons=P_sale_state_cons.*mult;


   [M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state, Skwh_state, Scost_state, P_state_cons, P_sale_state_cons, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NP"));
 writematrix(M_sim, title);
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NP")); % damages
 writematrix(D_sim, title); 
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NP")); % buildings
 writematrix(BI_sim, title); 
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NP")); % Fiscal cost
 writematrix(Exp_sub_sim, title);


% do prices then subsidies
net_metering=0; % assume all states offer net metering

mult=1; % scale by certain value to make total installations equal to baseline
position=1; %which object to scale? 

 f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

    [mult, fval, exitflag, output]...
        = fminunc(f,mult,options);
    
Sunit_state_cons=Sunit_state_cons.*mult;
Skwh_state_cons=Skwh_state_cons.*mult;
Scost_state_cons=Scost_state_cons.*mult; 


    [M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state_cons, Skwh_state_cons, Scost_state_cons, P_state_cons, P_sale_state_cons, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NPS"));
 writematrix(M_sim, title);
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NPS")); % damages
 writematrix(D_sim, title); 
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NPS")); % buildings
 writematrix(BI_sim, title);  
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NPS")); % Fiscal cost
 writematrix(Exp_sub_sim, title); 

%6 XXXX do sun only
net_metering=net_metering_o; %reset to what it was before for decomp

mult=1; % scale by certain value to make total installations equal to baseline
position=4; %which object to scale? 4 is sunlight
 
 f = @(mult_choice)calc_decomp_dist(mult_choice, position, total_installs_base,params, A_t_cons, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

[mult, fval, exitflag, output]...
    = fminunc(f,mult,options);
    
A_t_cons=A_t_cons.*mult;

    [M_sim,BI_sim, D_sim,  DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
        simulate(params, A_t_cons, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

title=(sprintf('Outputs/M_t%s%s.csv', model_type,  "NSun"));
 writematrix(M_sim, title);
title=(sprintf('Outputs/D%s%s.csv', model_type,  "NSun")); % damages
 writematrix(D_sim, title);
title=(sprintf('Outputs/BI%s%s.csv', model_type,  "NSun")); % buildings
 writematrix(BI_sim, title); 
title=(sprintf('Outputs/Fisc%s%s.csv', model_type,  "NSun")); % Fiscal cost
 writematrix(Exp_sub_sim, title);


net_metering=net_metering_o; %reset to what it was before for decomp
 
end




  
 %% Optimal Policy 
 
 if run_opt_policy==1
     
     
spec="Opt";

if new_store==1
 spec=strcat(spec, "NS");

 frac_store=.75

 spec=strcat(spec,string(100*frac_store));
end

if nat_level==1
    spec=strcat(spec, "Nat");
end

if tract_level==1
    spec=strcat(spec, "TL");
end

if cost_neut==1
    spec=strcat(spec, "CN");
end

if net_cost==1
    spec=strcat(spec, "NC");
end

if cost_half==1
    spec=strcat(spec, "CH");
end

if min_damages==1
    spec=strcat(spec, "MD");
end

if mcf150==1
    spec=strcat(spec, "mcf150");
end

if mcf125==1
    spec=strcat(spec, "mcf125");
end

if clean_coal==1
    title='../Data/CleanedData/ProdMedianCC.csv';
    ProdMedian=readmatrix(title);
    
   % title='../Data/CleanedData/DamagesCoef1CC.csv';
   
    title=(sprintf('../Data/CleanedData/DamagesCoef1CC%d.csv',scc)); 
    %end
    DamagesCoef1=readmatrix(title);
    
    % put into 1000s of dollars
    DamagesCoef1=DamagesCoef1/1000;
    
    %title='../Data/CleanedData/DamagesCoef2CC.csv';

   
    title=(sprintf('../Data/CleanedData/DamagesCoef2CC%d.csv',scc)); 
    
    DamagesCoef2=readmatrix(title);
    
    DamagesCoef2=DamagesCoef2/1000;
end


%increase_renew
   if increase_renew==1
        for r_i=1:nnerc
            title=(sprintf('../Data/CleanedData/RenProdL_r%d.csv',r_i));    
            RenProd(:, r_i)=readmatrix(title);   
        end
   end

   if increase_renew==2
        for r_i=1:nnerc
            title=(sprintf('../Data/CleanedData/RenProdM_r%d.csv',r_i));    
            RenProd(:, r_i)=readmatrix(title);   
        end
   end

    if increase_renew==3
        for r_i=1:nnerc
            title=(sprintf('../Data/CleanedData/RenProdH_r%d.csv',r_i));    
            RenProd(:, r_i)=readmatrix(title);   
        end
   end

     
 N_mean=15;
 run_damages=1;
 dW_ds_tol=.1;
 % dW_ds_tol=1

 if nat_level==1
     dW_ds_tol=dW_ds_tol./100;
 end

 if trans_con_loc==1
     dW_ds_tol=dW_ds_tol*3;
 end
     
     %run once to get baseline
     state_max=1;
   [M_sim,BI_sim, D_sim, DPol_sim, MargDamagesHour, M_bin_sim, pi_bin, Exp_sub_sim, Tot_ut_sim]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

   

     
    sub_amt_base_t=calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
    Sunit_t, Scost_t, Skwh_t, P_t, PSale_t, PIns_t, PIns0_t);  

    sub_amt_base_t=sub_amt_base_t(:,1); % just one per tract and all bins are identical here

    N_tract_state=sum((ones(1,nstate)).*tract_state_xxwalk_mat)'; 

    sub_base_state=sum((sub_amt_base_t*ones(1,nstate)).*tract_state_xxwalk_mat)';
    sub_base_state=sub_base_state./N_tract_state; %this one is an average, not a sum
    sub_base_nat=sum(sub_amt_base_t)./sum(N_tract_state);

    
    
    BI_base_state=sum((BI_sim*ones(1,nstate)).*tract_state_xxwalk_mat)'; 
    M_base_state=sum((M_sim*ones(1,nstate)).*tract_state_xxwalk_mat)'; 
    Fisc_base_state=sum((Exp_sub_sim*ones(1,nstate)).*tract_state_xxwalk_mat)'; 
   D_base_state=sum((D_sim*ones(1,nstate)).*tract_state_xxwalk_mat)';  % not clear ths

    D_base=D_sim;
    DPol_base=DPol_sim;
  %  Exp_sub_base=

     state_max=0 % subject to state level maximum subsidies from data
     

 % good guess    
 
     Sunit_guess=zeros(nstate,1); % inital guess
     Scost_guess=Sunit_guess;
     

     ResProd=calc_res_prod(A_t, M_sim, TractSunlightProfile); %yearly res prod in mwh

     marg_damages_t=calc_marg_damages4(MargDamagesHour, TractSunlightProfile, ResProd); % this is per year. multiplied by A_t later to get lifetime



 % cost in dollars of 1 unit increase in electricity in each tract
     
     HH_t=sum(HH_bin,2);
     
     ave_marg_damages_state=   sum((marg_damages_t.*HH_t*ones(1,nstate)).*tract_state_xxwalk_mat)' ... 
         ./ sum((HH_t*ones(1,nstate)).*tract_state_xxwalk_mat)';
     
     Skwh_guess=-ave_marg_damages_state; 

     if over_dist==1
         Skwh_guess=Skwh_state;
         Sunit_guess=Sunit_state;
         Scost_guess=Scost_state;
     end
 %}
  %bad guess   
 %   Skwh_guess=mean(ave_marg_damages_state)*ones(nstate,1); bad guess

 if marginal_increase==1 % change subsidies around current level
    Skwh_guess=Skwh_state;
    Sunit_guess=Sunit_state;
    Scost_guess=Scost_state;

    marginal_increase_global=1

 end

 if tract_level==1
     tl_global=1
    dW_ds_tol=50;


    %use state level optima for starting points for tract level
    if min_damages==0
        load(sprintf('Outputs/Data%s%s', model_type,  'OptCN'), 'Sunit_opt', 'Scost_opt', 'Skwh_opt')
    elseif min_damages==1
        load(sprintf('Outputs/Data%s%s', model_type,  'OptMD'), 'Sunit_opt', 'Scost_opt', 'Skwh_opt')
    end 

    Skwh_guess=Skwh_opt;
    Sunit_guess=Sunit_opt; %XXXX
    Scost_guess=Scost_opt;

 end



   RenProd_o=RenProd;
    

   if new_store==1


        RenProd_o=RenProd;

       RenTotYear=sum(RenProd,1);
       
       LoadTot=sum(LoadRaw,1);
       
       LoadProfile=LoadRaw./(ones(nhours,1)*LoadTot);
       
       RenProd=frac_store*(ones(nhours,1)*RenTotYear).*LoadProfile+(1-frac_store)*RenProd;
       
       for t_i=1:ntract
           
            reg=tract_nerc_xxwalk(t_i,1);
            TractSunlightProfile(t_i,:)=frac_store*LoadProfile(:,reg)'+(1-frac_store)*TractSunlightProfile(t_i,:);
       end
       
   end
   
   
   if cost_neut==0
       lambda=1

       if mcf150==1
          lambda=1.5
       end

       if mcf125==1
          lambda=1.25
       end
   end
   
   if cost_neut==1 || min_damages==1 || cost_half==1 || net_cost==1
      
       CostBase=sum(Exp_sub_sim);
             
       
       lambda_guess=.5433938;

       if scc==152 || scc==148
           lambda_guess=0.87;
       end

       if scc==330
           lambda_guess=1.3;
       end



        if nat_level==1
            
             lambda_guess=.5201;

           if scc==152 || scc==148
               lambda_guess=0.97;
           end

        end

        if clean_coal==1
            lambda_guess=0.503;
           if scc==152 || scc==148
               lambda_guess=0.97;
           end
        end


        if net_cost==1
            CostBase=CostBase+D_sim;
             lambda_guess=.44158;
            if scc==152 || scc==148
               lambda_guess=0.90;
           end
        end
       
       if cost_half==1
           CostBase=CostBase*.9;
     %      lambda_guess=.54568;
       end  
       
       if min_damages==1
           lambda_guess=.209;
            if scc==152 || scc==148
                lambda_guess=.43;
            end

            if scc==330
                lambda_guess=.75;
            end

       end


        if line_losses==1
            lambda_guess=.561;
            if scc==152 || scc==148
               lambda_guess=0.99;
           end
        end

        if nat_level==1 && min_damages==1
            
             lambda_guess=.157;
        end         

        if min_damages==1 && line_losses==1
            lambda_guess=.2;
        end

        if tl_global==1 && min_damages==0
            lambda_guess=.543; %XXX
            if scc==152 || scc==148
               lambda_guess=0.97;
           end
        end

        if tl_global==1 && min_damages==1
           % lambda_guess=.190;
            lambda_guess=.4358645599;
        end        

        if over_dist==1
            lambda_guess=0.872423046;

            if min_damages==1
                lambda_guess=0.43722541198;
            end
        end

        lambda_guess
       
 %      Cost_dist=  calc_cost_dist(CostBase, Sunit_guess,Scost_guess,Skwh_guess,  params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
 %P_state, P_sale_state, PIns_state,PIns0_state, nat_level, lambda_guess, min_damages, net_cost, dW_ds_tol, ...
 %   Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 %ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);

 
 %turn on to get fast initial guess with lower tolerance
 
 
     dW_ds_tol=dW_ds_tol*10; % looser for speed then re-tighten
     
 % get starting guess
 if tl_global==0
    [Sunit_guess, Skwh_guess, Scost_guess]=solve_S2(Sunit_guess,Scost_guess,Skwh_guess,  params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
 P_state, P_sale_state, PIns_state,PIns0_state, nat_level, lambda_guess, min_damages, net_cost, dW_ds_tol, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);
 end

options = optimset...
      ('Display','iter-detailed', 'MaxIter', 10, 'MaxFunEvals', 20, 'TolFun', 1E-3);
  

    
    f = @(lambda_choice)calc_cost_dist(CostBase, Sunit_guess,Scost_guess,Skwh_guess,  params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
 P_state, P_sale_state, PIns_state,PIns0_state, nat_level, lambda_choice, min_damages, net_cost, dW_ds_tol,...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);


    if over_dist==0 && (tl_global==0 || min_damages==0)
        [lambda_guess, fval, exitflag, output]...
            = fminsearch(f,lambda_guess,options);

    end
    
    dW_ds_tol=dW_ds_tol/10;
 
    
        f = @(lambda_choice)calc_cost_dist(CostBase, Sunit_guess,Scost_guess,Skwh_guess,  params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    P_state, P_sale_state, PIns_state,PIns0_state, nat_level, lambda_choice, min_damages, net_cost, dW_ds_tol,...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
    ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);


options = optimset...
      ('Display','iter-detailed', 'MaxIter', 20, 'MaxFunEvals', 30, 'TolFun', 1E-3);

    [lambda, fval, exitflag, output]...
        = fminsearch(f,lambda_guess,options);

%lambdaoff=1
%lambda=lambda_guess

   end 





if new_store==1 && cost_neut==0 && min_damages==0 && no_reopt==1 % just run with current subsides but new storage capacity
    Sunit_opt=Sunit_state;
    Skwh_opt=Skwh_state;
    Scost_opt=Scost_state;

    state_max=1
else
        [Sunit_opt, Skwh_opt, Scost_opt]=solve_S2(Sunit_guess,Scost_guess,Skwh_guess,  params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
     P_state, P_sale_state, PIns_state,PIns0_state, nat_level, lambda, min_damages, net_cost, dW_ds_tol, ...
        Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
     ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile);


end
  
    
        [M_opt,BI_opt, D_opt,  DPol_opt, MargDamagesHour_opt, M_bin_opt, pi_bin_opt, Exp_sub_opt, Tot_ut_opt]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_opt, Skwh_opt, Scost_opt, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);
 
 
%once to get damages offset (so damages with no rooftop solar and no
%technology improvements
[D0, DPol0, ~]=calc_damages5(0*ResProd, Load, RenProd_o, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ProdMedian, ...
    DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2); % damages\\

ResProd=calc_res_prod(A_t, M_opt,TractSunlightProfile); % this uses data level installations. Could also do with predicted installations


Eload=Load-ResProd-RenProd;

exceedEload=(Eload>ELoadMax);

meanExceed=mean(exceedEload,'all')

exceedEload99=(Eload>ELoadMax99);

meanExceed99=mean(exceedEload99,'all')

 if no_reopt==1
  model_type=strcat(model_type, "NR");
 end


 if line_losses==1
  model_type=strcat(model_type, "LL");
 end

  if trans_con==1
  model_type=strcat(model_type, "TC");
 end

if increase_renew==1
    model_type=strcat(model_type, "IRL");
end

if increase_renew==2
    model_type=strcat(model_type, "IRM");
end

if increase_renew==3
    model_type=strcat(model_type, "IRH");
end


if clean_coal==1
    model_type=strcat(model_type, "CC");
%    spec=strcat(spec, "CC");
end





if scc~=148
    spec=sprintf('%s%d', spec, scc);
%    spec=strcat(spec, "CC");
end

lambda

% model_type=strcat(model_type,string(scc));

clear  TractSunlightProfile PPRandStacked StateSunlightProfile Nbar_i    M_bin_opt M_bin_sim ...
    A_t M_t BI_t Nbar_i HH_bin CollPerc_t Pol_t Own_t tract_nerc_xxwalk tract_state_xxwalk pi_bin pi_bin_opt ...
    statemaxtaxcredit_t statetaxcredit_t sunit_cap_t Scost_t Skwh_t Sunit_t PPRegion PPRegion2 Select MargDamagesHour ...
    MargDamagesHour_opt RenProd ResProd Load LoadRaw 

save(sprintf('Outputs/Data%s%s', model_type,  spec), 'Sunit_opt', 'Scost_opt', 'Skwh_opt', 'N_mean', 'ntract', ...
    'P_t', 'PSale_t', 'PIns_t', 'PIns0_t', 'nstate', 'N_tract_state', ...
'BI_opt', 'M_opt', 'Exp_sub_opt', 'D_opt', 'DPol_opt', 'HH_t', 'sub_base_state', 'BI_base_state', 'Fisc_base_state', 'D_base_state', ...
'sub_base_nat', 'BI_base_state', 'Fisc_base_state','D_base', 'DPol_base', 'D_base_state', ...
'M_sim','BI_sim', 'D_sim', 'DPol_sim', 'Exp_sub_sim', 'D0', 'DPol0', 'Tot_ut_sim' , 'Tot_ut_opt', 'meanExceed', 'meanExceed99');


if tl_global==1
    save(sprintf('Outputs/TractLevel%s%s', model_type,  spec), 'Sunit_tg', 'Skwh_tg', 'Scost_tg');
end


 end
 
 
 
 done=1